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Abstract 

This paper contains an analysis of multimode solutions expressible in terms of Riemann invariants 
obtained from complex characteristic elements. A variant of the conditional symmetry method for 
constructing this type of solution is proposed. Examples of applications of the proposed approach 
to a nonstationary irrotational ideal plastic flow in (2+1) dimensions are discussed in detail. Several 
new classes of solutions, some of them bounded, have been obtained in closed form. For a particular 
case of a stationary ideal plastic flow, extrusion dies have been drawn. 
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1 Introduction 

Riemann waves represent a very important class of solution of nonlinear first-order systems of PDEs. 
They are ubiquitous in the differential equations of mathematical physics, since they appear in all multi- 
dimensional hyperbolic systems and constitute their elementary solutions. Their characteristic feature is 
that, in most cases, they are expressible only in implicit form. For a homogeneous hyperbolic quasilinear 
system of first-order PDEs, 

du a 

A %( u )~faT =0 > i = 1 > ••->£> a = l,...,q, fjL=l,...,m, (1) 

(where A 1 , . . . , A p ) are q x m matrix functions of an unknown function u. Here we adopt the convention 
that repeated indices are summed unless one of them is in a bracket. 

A Riemann wave solution is defined by the equation u = f(r(x, u)), where / : R — > W 1 and the function 
r(x, u) = Xi(u)x l is called the Riemann invariant of the vector A satisfying the equation ker (A^A^u)) ^ 
0. These solutions have rank at most equal to one. They are building blocks for constructing more 
general types of solutions describing nonlinear superpositions of many waves (k- waves), which are very 
interesting from the physical point of view. Until now, the only way to approach this task was through 
the generalized method of characteristics (GMC) [3j [7l ETJ |22j [29J, [30] and more recently through the 
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conditional symmetry method (CSM) [U0H31HI3IIE1I2E]- The GMC relies on treating Riemann invariants 
as new dependent variables (which remain constant along the appropriate characteristic curves of the 
initial system |lj and constitute a set of invariants of the Abelian algebra of some vector field X a = (it)<9* 
with Xf = for 1 < a < k < p. This leads to the reduction of the dimension of the problem. The 
most important theoretical results obtained with the use of the GMC or CSM [T^ include the finding 
of necessary and sufficient conditions for the existence of Riemann k- waves in multidimensional systems. 
It was shown |29j that these solutions depend on k arbitrary functions of one variable. Some criteria 
were also found |29| for determining the elastic or nonelastic character of the superposition between 
Riemann waves described by hyperbolic systems, which is particularly useful in physical applications. 
In applications to fluid dynamics and nonlinear field theory, many new and interesting results were 
obtained [HI HOI HU H21 12U HH1 IM H21 EH • Both the GMC and CSM methods, like all other techniques 
for solving PDEs, have their limitations. This fact has motivated the authors to search for the means 
of constructing larger classes of multiple wave solutions expressible in terms of Riemann invariants by 
allowing the introduction of complex integral elements in the place of real simple integral elements (with 
which the solutions of hyperbolic systems are built [29 ). This idea originated from the work of S. 
Sobolev |31| in which he solved the wave equation by using the associated complex wave vectors. We are 
particularly interested in the construction of nonlinear superpositions of elementary rank— A: solutions 
(modes) and the proposed analysis indicates that the language of conditional symmetries is an effective 
tool for this purpose. This approach is applied to the nonstationary irrotational flow of a ideal plastic 
material in its elliptic region. 

The organization of this paper is as follows. Section 2 contains a detailed account of the generalized 
method of characteristics for first-order quasilinear systems of PDEs in many dimensions based on complex 
characteristic elements. In section 3 we formulate the problem of multimode solutions expressible in terms 
of Riemann invariants by means of a group theoretical approach. This allows us to formulate the necessary 
and sufficient conditions for constructing these types of solutions. Section 4 presents an application of 
the developed approach to the system of equations describing a nonstationary irrotational flow of ideal 
plastic materials in (2+1) dimensions. Several new classes of solutions are obtained in closed form. We 
obtain several bounded solutions. For the stationary case, limiting ourselves to the region where the 
gradient catastrophe does not occur, we have drawn extrusion dies and the flow inside the tools. Section 
5 summarizes the obtained results and contains some suggestions regarding some further developments. 

2 Complex characteristic elements. 

The methodological approach assumed in this section is based on the generalized method of characteristics 
which has been extensively developed (e.g. in [3, 9lll4 [ [T7 j l29| and references therein) for multidimensional 
homogeneous and nonhomogeneous hyperbolic systems of first-order PDEs. A specific feature of that 
approach is an algebraic and geometric point of view. An algebraization of systems of PDEs was made 
possible by representing the general integral elements as linear combinations of some special elements 
associated with those vector fields which generate characteristic curves in the spaces of independent 
variables X and dependent variables U, respectively. The introduction of these elements (called simple 
integral elements) proved to be very useful for constructing certain classes of rank-fc solutions in closed 
form. These integral elements proved to correspond to Riemann wave solutions in the case of nonclliptic 
systems and serve to construct multiple waves (A:-waves) as a superposition of several single Riemann 
waves. 
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The generalized method of characteristics for solving quasilinear hyperbolic first-order systems of can 
be extended to the case of complex characteristic elements |29| . These elements were introduced not 
only for elliptic systems but also for hyperbolic systems by allowing the wave vectors to be the complex 
solutions of the dispersion relation associated with the initial system of equations. 



The starting point is to make an algebraization, according to 0E, 29 , of a first-order system of PDEs 
in p independent and q dependent variables written in its matrix form 

A l {u)ui = Q, i = l,...,p, 

x = (x 1 ,...,x p ) e X C RP, u= (u\...,u q ) eUCRQ, 

where A A p are to x q matrix functions of u and we denote the partial derivatives by uf = du a jdx 1 . 
The matrix Lf satisfying the conditions 

ufe{Lf:A i *Lf = 0, s = l,...,q} 

at some open given point u$ G U is called an integral element of the system This matrix L = 
[du a / dx 1 ) is a matrix of the tangent mapping du : X — > T U U given by the formula 

X 3 (5x l ) -> (Su a ) e T U U, where Su a = ufSx 1 . 

The tangent mapping du(x) determines an element of linear space L(X,T U U), which can be identified 
with the tensor product T U U <g> X* , (where X* is the dual space of X, i.e. the space of linear forms). It 
is well known [31 Q31 that each element of this tensor product can be represented as a finite sum of 
simple tensors of the form 

L = 7® A, 

where A € X* is a covector and 7 = 7 a gf^ £ T U C/ is a tangent vector at the point u G U. Hence, the 
integral element Lf is called a simple element if rank(L") = 1. To determine a simple integral element 
Lf we have to find a vector field 7 g T U U and a covector A € X* satisfying the so-called wave relation 

(X i A i (u))r = 0- (3) 

The necessary and sufficient condition for the existence of a nonzero solution 7 for the equation (|3| is 

rank (\iA l (u)) < min(m, q). (4) 

This relation is known as the dispersion relation. If the covector A = \i(u)dx l satisfies the dispersion 
relation Q then there exists a polarization vector 7 € T U U satisfying the wave relation (|3j. The algebraic 
approach which has been used in [3J Q31 [53] for hyperbolic systems of equations ([2| allows the construction 
of certain classes of fc-wave solutions admitting k arbitrary functions of one variable. The replacement 
of the simple real element Lf by the matrix of derivatives uf in the system of equations Q allows us to 
construct more general classes of solutions by replacing the real elements with complex ones. A specific 
form of solution u(x) of ^ is postulated for which the tangent mapping du(x) is a sum of a complex 
element and its complex conjugate 

du a (x) = £(x)7 Q (w)A 4 (w)da; 1 + i(x)f(u)X i (u)dx i , (5) 
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where 7 = (7 1 , ...j'JeC and A = (Ai, . . . , A p ) G C p satisfy 

A^ S (")7 Q = 0, \iA™(u)j a =0. 

Here the quantity £(x) 7^ is treated as a complex function of the real variables x. In what follows we 
assume that the vectors 7 and 7 are linearly independent. The proposed form of solution |5]) is more 
general than the one proposed in [18 , 29 for which the derivatives uf are represented by a real simple 
element leading to a simple Riemann wave solution. To distinguish this situation with the one proposed 
in ([5]), we call the solution associated with a complex element and its complex conjugate a simple mode 
solution. Similarly, as in the case of /c-waves we have to find the necessary and sufficient conditions 
for the existence of solutions of type (pjj. These conditions are called involutivity conditions. 



First we derive a number of necessary conditions on the vector fields (7, A) and their complex conjugate 
(7, A) as a requirement for the existence of rank-2 (simple mode) solutions of the original system 
Namely, closing (|5| by exterior differentiation, we obtain the following 2-forms 

j a (d£ A A + £d\) + fd 7 A A + 7 Q (dC A A + £dX) + fjfy A A = (6) 

which have to satisfy Using ([5| we get 

rfA = (A A A :7 + £A A A i7 , dj = £7^ <g> A + 7^ ® dX, (7) 

where we have used the following notation 

= 7 a^ A ' 7 - = 7 s^ 7 - 

Substituting ^ into the prolonged system ^ we obtain 

7 [d( A A + ££A A A )7 ] + 7 <g> [d£ A A + ££A A A, 7 ] + [7, 7] <g) A A A = 0, (8) 
whenever the differential §5§ holds. The commutator of vector fields 7 and 7, is denoted by 

[7, 7] = (7, l)u + 7,AiAi, 7 - 7,Ai Ai, 7 , 

while by (7,7),. we denote a part of the commutator which contains the differentiation with respect to 
the variables u a , i.e. 

d d 
(7,7)„ =7°^ 7(«,A)| A=const . - 7 Q — 7(u,A)| x=conilt . . 

Let $ be an annihilator of the vectors 7 and 7, i.e. 

< ujjj >— 0, < wj7 >= 0, we$ = An{7,7}. 

Here, by the parenthesis < luj^ >, we denote the contraction of the 1-form uj £ T*U with the vector field 
7 € T„t7 and similarly for the vector field 7 <G T u {7. Multiplying the equations Q by the 1-form lj € $ 
we get 

< wj [7, 7] > A A A = 0. (9) 
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We look for the compatibility condition for which the system (|9J) does not provide any algebraic constraints 
on the coefficients £ and £. This postulate means that the profile of the simple modes associated with 
7 (8 A and 7 ® A can be chosen in an arbitrary way for the initial (or boundary) conditions. It follows 
from ([9]) that the commutator of the vector fields 7 and 7 is a linear combination of 7 and 7, where the 
coefficients are not necessarily constant, i.e. 

[7,7] £ span {7, 7} . 

So, the Frobenius theorem is satisfied. That is, at every point uo of the space of dependent variables U , 
there exists a tangent surface S spanned by the vector fields 7 and 7 passing through the point uo £ U . 
Moreover the above condition implies that there exists a complex- valued function a such that 

[7, 7] =aj -«7, (10) 

since 

[7)7] = [7>7] = - [7)7] 

holds. 



Next, from the vector fields 7 and 7 defined on [/-space, we can construct the coframe <J/, that is the 
set of 1-forms <j\ and a-i defined on U satisfying the conditions 



< o 1 J7 >= 1, 

< <t 2 j7 >= 0, 



< <j 1 J7 >= 0, 

< CT 2 J7 >= 1. 



0-1,0-2 e 



(11) 



Substituting ^ and ( 10 1 into the prolonged system (JsJ) and multiplying by o~\ and 02 respectively, we 

(i) <i£ A A + £ (£A A A )7 + £A A A )7 ) + a££A A A = 0, 
(m) c^AA + |(CAA A 7 + fAAA, 7 ) + q^AA = 0. 



If d£, = then equation (12 i) is a complex conjugate of (|12|ii). So one can consider one of them, 
say (12 i). We look for the condition of integrability such that the system (12 i) does not impose any 
restriction on the form of the coefficients £ and £. By exterior multiplication of ( 12 i) by A, we obtain 



A 7 A A A A 







and its respective complex conjugate equation 



A, 7 A A A A = 0. 



(13) 



(14) 



Note that if conditions (10 1, ([13]) and ( 14 1 are satisfied, then from the Cartan lemma the system ( 12 1 and 
consequently the set of 1-forms Q has nontrivial solutions for d£ and So under these circumstances 
the following result holds |29| . 

Proposition 1 The necessary condition for the Pfaffian system |5p to possess a simple mode solution is 
that the following two constraints on the complex-valued vector fields 7 and A be satisfied 



(i) [7, 7] = Q-7 - o<7, 



for any complex-valued function a£C 



(m) A 7 , 7 j7 £ span {A, A} 
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Note that the conditions ( 10 1 , ( 13 ) and ( 14 1 are strong requirements on the functions 7 = (7 1 , . . . , 7 9 ) 
and A = (Ai, . . . ,X p ) and their complex conjugates which satisfy the wave relation Consequently, 
the postulated form of a one-mode solution u(x) of the elliptic system ([2j, required by the generalized 
method of characteristics (GMC) |21[ [29] . is such that all first-order derivatives of u(x) with respect 
to x l are decomposable in the form ([5]). This restriction is a strong limitation on the admissible class 
of solutions of pj. So far, there have been few, if any, known examples of such solutions. Therefore, 
it is worth developing this idea by weakening the integrability conditions (10 1, (13 1 and (14 1 in order 
to construct multi-mode solutions. In the next section we propose an alternative way of constructing 
multi-mode solutions (expressed in terms of Riemann invariants) which are obtained from a version of 
the conditional symmetry method (CSM) |16| by adapting it to the elliptic systems. This constitutes the 
objective of this paper. 



3 Conditional symmetry method and multimode solutions in terms 
of Riemann Invariants. 

In this section, we examine certain aspects of the conditional symmetry method in the context of Riemann 
invariants for which the wave vectors are complex solutions of the dispersion relation Q associated with 
the original system ^ . Let us consider a nondegenerate first-order quasilinear system of PDEs (J2J) in its 
matrix form 

A 1 (u)u 1 + ... + A p (u)u p = 0, (15) 

where A 1 , . . . , A p are m x q real- valued matrix functions of u. Let us suppose that there exist k linearly 
independent non-conjugate complex-valued wave vectors 

\ A (u) = (Xf{u),...,X A (u)) 6P, X A ^X B V B e{l,...,k}, A = l,...,k <p (16) 



which satisfy the dispersion relation (|4j. One should note that in (16 1 we do not require indices A=^ B, 
which means that real wave vectors are excluded from our consideration (we do not consider here the 
mixed case for wave vectors involving real and complex wave vectors). Under the above hypotheses, k 



wave vectors ( 16 ) and their complex conjugates 

X A (u) = (X A («), . . . , X A (u)) 6P, A = l,...,k, 

satisfy the dispersion relation Q. In what follows it is useful to introduce the notation c.c. which means 
the complex conjugate of the previous term or equation. This notation is convenient for computational 
purposes allowing the presentation of some expressions in abbreviated form. 

Let us suppose that there exists a unique solution u(x) of the equations ([TJ of the form 

u = f(r 1 (x, u), . . . , r k (x, u), f 1 (a;, u), . . . , f k (x, it)) + c.c, (17) 

where the complex- valued functions r A ,r A : W x E 9 — > C are called the Riemann invariants associated 
respectively to wave vectors X A , X A and are defined by 

r A {x,u) = X A {u)x\ f A (x,u) = X A (u)x\ A = l,...,k. (18) 
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Note that the functions u(x) are defined implicitly in te rms of u a ,x t , r A and r A . For any function 



/ : C k — > C 9 and its complex conjugate, the equation (17) determines a unique real function u(x) on a 
neighborhood of the origin x = 0. The Jacobi matrix of derivatives of u(x) is given by 



df df\ dr 



dr dr I du 



du = «) = (J, - 
or equivalently as 

1 dr ' dr J 

where the k x k matrices M 1 and M 2 are defined by 



+ c.c. 



df + df 
dr dr 



A + c.c. 



(19) 



(20) 



M 1 



dr f df df 



du \ dr dr 

d L /dl + df 
du \ df df 



h- 



M 
dr 



du \ dr dr ) ^ 
\ dr A J dr 



df/dl + df 
du \ df df 

df ( df df 



df (df 
du \ dr 
-l 



df 
dr 



du \ dr dr 

eC qxk , A = (xf) e C kxp , 



(21) 



or 

df A 



dr 
du 



(22) 



(23) 



and their respective conjugate equations. The matrices I q and Ik are the q x q and k x k identity matrices 
respectively. We use the implicit function theorem to obtain the following conditions ensuring that r A 7 f A 
and u° are expressible as graphs over some open subset D C K p , 



or 



det ( Jfc 



det ( Ig - 



dr ( df df 



df df \ dr 



dr dr I du 



+ c.c. 



du V dr dr 



dr(dl + dl 
du \ dr dr 



h- 



7^ 0, and det 



dr ( df df 



^0 



dr(dl + df 
du \ df df 



du V dr dr 



5r (df_ df_ 
du \ df df 



(24) 



(25) 



7^0 



So the inverse matrix in ( 19 1 or in (20 1 is well-defined in the vicinity of x = 0, since 

dr 



du 



= o, 



dr 

du 



= at x = 0. 



(26) 



Note that on the hypersurface defined by equation (17) and one of the expressions (24 1 or (25 1 equal to 



zero, the gradient of the function u{x) becomes infinite for some value of x % . So the multimode solution, 
given by expression (17) loses its sense on this hypersurface. Consequently, some types of discontinuities, 



i.e. shock waves can occur. In what follows, we search for solutions defined on a neighborhood of x = 



7 



under the assumption that the conditions (24 1 or ( 25 1 are different from zero. This means that if the 



initial data is sufficiently small, then there exists a time interval [io,? 1 ], T > to (where we denote the 
independent variables by t = 71 where p = n + 1) in which the gradient catastrophe for 

the solution u(t,x) of the system ([2| does not take place |30| . 

Note that the Jacobian matrix of u(x) has at most rank equal to 2k. It follows that the proposed 



solution (17 1 is also at most of rank-2/c. Its image is a 2fc-dimensional submanifold §2k m the first jet 
space J 1 =J 1 (XxU). 



Let us introduce now a set of p — 2k linearly independent vectors £ a 



defined by 



&(«) = (£(«),... 

satisfying the orthogonality conditions 



,g(u)), a=l,...,p-2k, 



A=l, 



. , k, a = 
i A fc , A 1 , . . 



l,...,p-2k, 



(27) 



(28) 



, A fc } . It should be noted that the 



for a set of 2k linearly independent wave vectors {A 1 , 
set {A 1 , . . . , A fc , A 1 , . . . , A fe , . . . , £P~ 2k } forms a basis for the space of independent variables X. Note 
also, that the vectors £ a are not uniquely defined since they obey the homogeneous conditions ( 28 ) . As a 



consequence of equation (19 I or (20 1, the graph T = {x,u(x)} is invariant under the family of first-order 
differential operators 

X = £(u)^, a = l,...,p-2k, (29) 

defined on X x U space. Since the vector fields X a do not include vectors tangent to the direction of u, 
they form an Abelian distribution on X x U space, i.e. 



[X a ,X b ] = 0, a,b=l,...,p-2k. 



(30) 



The set {r 1 , . . . , r k , r 1 , . . . , r k , u 1 , . . . it 9 } constitutes a complete set of invariants of the Abelian algebra 
C generated by the vector fields (29 1. So geometrically, the characterization of the proposed solution (17 1 



of equations ( |15| can be interpreted in the following way. If u(x) is a g-component function defined on a 
neighborhood of the origin x — such that the graph of the solution T = {(x, u(x))} is invariant under a 



set of p — 2k vector fields X a with the orthogonality property (28 1, then for some function / the expression 



u(x) is a solution of equation (17 1. Hence the group-invariant solutions of the system ( 15 1 consist of those 
functions u — f(x) which satisfy the overdetermined system composed of the initial system (15) together 
with the invariance conditions 



£>? = 0, i = l,...,p, a = l,...,p-2k, 
ensuring that the characteristics of the vector fields X a are equal to zero. 



(31) 



It should be noted that, in general, the conditions (311 are weaker than the differential constraints 
^ required by the generalized method of characteristics, since the latter method is subjected to the 
algebraic conditions |3|. In fact, equations (31 1 imply that there exist complex- valued matrix functions 
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&a(x, u) and & A (x, u) defined on the first jet space J 
respect to x % are decomposable in the following way 



J(X x U) such that all first derivatives of u with 



u? = $ a A {x,u)\ A +W* A (x,u)\i 



(32) 



where 



n=[r q - 



9/ <9/\ dr 



<9r 3r / 9u 
df df\ dr 



dr dr J du 



+ c.c. 



9f + df 
dr dr 

df df 



(33) 



n = 



jr + ir) M 

dr dr J 

df df\ —1 

+ IF M 
or ar / 



df df \ —2 

IF + IF I M 
or or 

df df \ , , 
or or 



(34) 



The matrices ^^A^ and ^^Xf appearing in equation (32 1 do not necessarily satisfy the wave relation 
|3]). As a result of this fact, the restrictions on the initial data at t = are eased, so we are able to 
consider more diverse types of modes in the superpositions than in the case of the generalized method of 
characteristics described in section 2. 



constraints (31 ) 



Let us now proceed to solve the overdetermined system composed of equations (15 1 and differential 

(35) 



Ag(u)u? = 0, C a (u)uf = 0. 



Substituting (19 1 or (20 1 into (15 1 yields the trace condition 

-1 



tr UM I q - 



df df\ dr 



dr dr I du 



+ c.c. 



dl + df 
dr dr 



or 



tr A" 



dr 



+ ~cfr) ( M X + M X )+ c - c 



X + c.c. 



= 



= 



(36) 



(37) 



on the wave vectors A and A and on the functions / and /, where A , . . . , A q are p x q matrix functions of 
u (i.e. A^ — (A^(u)) £ M. pxq , \i = 1, . . . , m). For the given initial system of equations ( 15 1, the matrices 



A 11 are known functions of u and the trace conditions (36 1 or (37 1 are conditions on the functions /, /, A 



A (or on £ due to the orthogonality conditions (28 1). From the computational point of view, it is useful 
to split x l into x lA and x la and to choose a basis for the wave vector X A and X A such that 



A = dx l 



X? dx l 



X = dx l 



Xfdx 1 



A = l, 



, fc, 



where (iA,i a ) is a permutation of (1, . . . ,p). So, the expressions (23 1 become 

dX A , 



dr A 
du a 



du a 



df A 
du a 



du a 



(38) 



(39) 
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Substituting (39 I into (36) (or (37 1) yields 



tr (l q - Q a x % « 



or 



dR dR 



1,1 '"''l + I)^-^' 



A 



0, fi = 1 . . . , m, 



A = 0, ii = 1 



,m, 



(40) 
(41) 



where = (r , . 



h -1 



r ) and 



Qa = 

K a = 



df 


v 8f\ 




i%- 


_m 






df 


dr 


dr J 


du 




dr J 


du 




OR 



d*j a (df , df 

du \ dr Or 

dha (df , df 

du \ dr dr 



d*i a (dl + d£ 

du \ dr dr 

d^ a (dl.dp 

du \ dr dr . 



Va{ 0R + 8R ] 



(42) 



and for simplicity of notation, we note 





e C 2kxq , 



dR 



dj_ dj_ 
dr ' df 



G C 9X2fc , 



(43) 



for %a fixed and i a = 1, . . . ,p— 1. In (42 1 the 2k x 2k matrix K a is defined in terms of the k x k subblocks 
of the form d\i a /du (df /dr + d f / dr)/ ^ where r\ a is a matrix form of the block d\ a /du over the block 
dXa/du. The notation (df /dr, df /df) represents the matrix formed of the left block df /dr and the right 



block df /dr. Note that the functions 



and x %a are functionally independent in a neighborhood of 



x = and the matrix functions A 11 , df /dr, df /dr, df /dr, df /dr, Q a and K a depend on r and r only. 
So, equations (40) (or (41 1) have to be satisfied for any value of coordinates x ta . As a consequence, we 



have some constraints on these matrix functions. From the Cayley-Hamilton theorem, we know that for 
any n x n invertible matrix M, the expression (M~ l det M) is a polynomial in M of order (n— 1). Thus, 
using the tracelessness of the expression A 11 (l q — Q a x la ) (df/dR)A, we can replace equations (40 1 by 
the following condition 



0. 



where Q = adj(7 g — Q a x la ) € 



^qxq 



(44) 



where adj M denotes the adjoint of the matrix M. As a consequence the matrix Q is a polynomial of 
order (q — 1) in x' a . Taking (44 1 and all its partial derivatives with respect to x %a (with r, f fixed at 
x = 0), we obtain the following conditions for the matrix functions f(r, f) and X(f(r,f)) 



tr U"0( Ol • ■ • Qa s 



dl 
dR 



dl 
dR 



A 



0. 



(45) 



(46) 
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where s = 1, . . . , q — 1 and (ai, . . . , a s ) denotes the symmetrization over all indices in the bracket. A 
similar procedure can be applied to system (411 to yield (45 1 and 



tr A" 



d L + d l 
dR dR 



(ai 



Ka s )A) =0, 



(47) 



where now s = 1, 



2k — 1. Equations (45 1 represent the initial value conditions on a surface in the 

0. Note that equations (46 1 (or ( 47 1 ) form the conditions 



space of independent variables X, given at x 1 
required for the preservation of the property ( 45 1 along the flows represented by the vector fields ( 29 ) 

Xfd, 



Equation ( 38 1 allows us to express X a in the form 



X„ 



9r a 



A = l, 



,k. 



Substituting expressions (|42| into (|46| or (|47| and simplifying gives the unified form 
tr ( A" 



df df 
dR dR 



df df 
dR dR 



■la,) 



df df 

dR dR 



A 



0, 



(48) 



(49) 



The index s is either max(s) = q— 1 or max(s) = 2k — 1, we choose the one which is more convenient from 
the computational point of view. In this case, for k > 1 the two approaches, CSM and GMC, become 
essentially different and, as we demonstrate in the following examples, the CSM can provide rank-2fc 
solutions which are not Ricmann 2/c-waves as defined by the GMC since we weaken the integrability 
conditions Q for the wave vectors A and A. 

A change of variable on X x U allows us to rectify the vector fields X a and considerably simplify the 
structure of the overdetermined system ( 35 1 which classifies the preceding construction of multimodes 
solutions. For this system, in the new coordinates, we derive the necessary and sufficient conditions 
for the existence of rank-2fc solutions of the form (17 1. Suppose that there exists an invertible 2k x 2k 
subblock matrix 

L = (Ai), l<A,B<2k, (50) 
of the larger matrix A 6 (J2 2kx P j then the independent vector fields X a can be written as 

X a = d x a+2 k - (L~ 1 )gA a 3 +2k d x A, 



(51) 



introduce new coordinate functions 



which have the required form (29) and for which the orthogonality conditions (28 1 are fulfilled. We 

(52) 



z 1 = r : (x, u), . . . , z k — r k (x, u), z k+1 = f l {x, u), . . . , z 2k = r k {x, u), 



z 



2k+l 



„2k+l 



on X x U space which allow us to rectify the vector fields (51 1. As a result, we get 

d 



X, = 



Q z 2k+l ' ' 

The p-dimensional submanifold invariant under X\ 



X p -2ki is defined by equations of the form 



z k rz\ 



,z k ) + c.c. 



(54) 
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for an arbitrary function / : X — > U. The expression (54) is the general solution of the invariance 
conditions 

v z 2k+i , . . . , v z p = 0. (55) 



In general, the initial system (15 1 described in the new coordinates (i,n)elx[/isa nonlinear system 
of first-order PDEs. of the form 



A i lt() dPdv^_ Q i- k + 1 

.■in, -.dv 13 



i = 2k + 1, 

We obtain the following Jacobi matrix in the coordinates (z, z, v) 



,2k, 
■ ,P- 



(56) 



s dvPdz* ' 



(57) 



whenever the invariance conditions (55 1 are satisfied. Appending to the system (56 1 the invariance 
condition (551, we obtain the quasilinear reduced system of PDEs 



tr ^"(«)(/„- 
dv 



dv dR\ 1 dv 



q dz dv ) dz 
dv 



A =0, 



fi = l,...,m, 



(58) 



dz 2k+l QzP 



= 



or 



tr(A»g(/^ 



(9u 9z 



A = 0, n = l, 



dv 



d z 2k+l ■ 



dv 
dzP 



(59) 



We now provide some basic definitions that we require for the use of the conditional symmetry method 
in order to encompass the use of Ricmann invariants. 



A vector field X a is called a conditional symmetry of the original system (15) if X a is tangent to the 
manifold § = §a H Sq, i.e. 

e r (sSiuW) s, (60) 



pr«* Q 

where the first prolongation of X a is given by 



pr wx a = x a -e a , u ^u?—, 

and the submanifolds of the solution spaces are given by 



SA = 



a = 1, . . . ,p — 2k 



(j, = 1, 



•}■ 



(61) 



(62) 



12 



and 



§ Q = {(x, U W) : £»< = 0, a = 1, . . . , q, a = 1, . . . ,p - 2k} 
As a consequence, we can define an Abelian Lie algebra. 



(63) 



An Abelian Lie algebra C generated by the vector fields X\, . . . , X v -ik is called a conditional symmetry 
algebra of the original system (15 1 if the conditions 



pr 



(1) X a (A l Ul ) 



0, a = l,...,p-2k, 



are satisfied. 



(64) 



Supposing that £, spanned by the vector fields X±, . . . , X p ^2k, is a conditional symmetry algebra of 



the system ( 15 1, a solution u — f(x) is said to be a conditionally invariant solution of the system ( 15 1 if 
the graph L = {(x, f(x))} is invariant under the vector fields Xi, . . . , A p _2fc- 



Proposition. A nondegenerate quasilinear hyperbolic system of first-order PDEs ( 15) inp independent 
variables and q dependent variables admits a p— 2k -dimensional symmetry algebra C if and only if (p—2k) 
linearly independent vector fields X\, . . . ,X p _2k satisfy the conditions (4-5) and (4-9) on some neighborhood 



of (xo,Uq) of S . The solutions of (15) which are invariant under the Lie algebra C are precisely k-mode 
solutions of the form (17). 



proof The proof of this proposition is essentially similar to that of the proposition in |16J . For the 
sake of completeness we now give it. Let us express the vector fields X a in the new coordinates (z, v) on 
X x U. Equations (53 1 and (61 1 imply that 



pr 



^X a = X a , a = l,...,p-2k. 



(65) 



The symmetry criterion that has to be satisfied for G to be the symmetry group of the overdetermined 
system (58 1 (or ( 59 1 ) requires that the vector fields X a of G satisfy X a (A) = 0, whenever equation (58 1 



(or ( 59 1 ) is satisfied. Thus the symmetry criterion applied to the invariance conditions (55 1 vanishes 



identically. Applying this criterion to the system (56 1 in the new coordinates, carrying out the differenti- 



ation and taking into account the conditions (45) and (49 1, we obtain the equations which are identically 
satisfied. 



The converse is also true. 



The assumption that the system (15 1 is nondegenerate means, according 

Therefore, 



to |27| . that it is locally solvable and is of maximal rank at every point (xq,Uq) € S. Therefore, the 
infinitesimal symmetry condition is a necessary and sufficient condition for the existence of the symmetry 



group G of the overdetermined system ( 35 1 . Since the vector fields X a form an Abelian distribution on 



,r fc ,r\ 



p-2k 



The invariants 



X x U, it follows that conditions (45) and (49 1 are satisfied. The solutions of the overdetermined system 
( 35 1 are invariant under the algebra C generated by the p—2k vector fields X\ , . . . , X 
of the group G of such vector fields are provided by the functions {r 1 . 
the general multimode solution of (15 1 takes the required form (|17|). 



-k 1 



'}. So 



□ 
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4 The ideal plastic flow. 



In this section we would like to illustrate the proposed approach for constructing multimode solutions 
with the example of the ideal nonstationary irrotational planar ideal plastic flow subjected to an external 
force due to a potential V. Under the above assumptions the examined model is governed by a quasilinear 
homogenous system of five equations in (2+1) dimensions of the form [5, 20, 24 



a) 



(9 X cos 29 + 6y sin 29) + p (V x — u t — uu x — vu y ) = 0, 



b) <7 y — (9 X sin 29 — 9 y cos 29) + p (V y — Vt — uv x 

c) (u y + v x ) sin 29 + (u x — v y ) cos 29 = 0, 

d) u x + v y = 0, 

e) u y -v x = 0. 



vv y ) = 0, 



(66) 



Equations (66 a) and (66 b) involve the independent variables a,9,u,v and the potential V is a given 
function of (t, x, y) € W\ The stress tensor is defined by the mean pressure a and 9 is the angle relative 
to the x-axis minus ir/4. Equation (66 c) represents the Saint- Venant- Von Mises plasticity equation. 
Equations ( 66 d) and ( 66 e) for the velocities u and v (along the cc-axis and y-axis respectively) correspond 
to the incompressibility and the irrotationality of the flow of the plastic material. The independent 
variables are denoted by (x l ) = (t,x,y) E X C M 3 and the unknown functions by (u a ) = (a,9,u,v) g 
t/cl 4 . The first-order partial derivatives of a with respect to the independent variables x and y are 
denoted a x and o y respectively. The first-order partial derivatives of u, v and 9 are denoted similarly. 



Before using the approach described in section [3] to obtain solutions of the system (66), we begin by 
expressing the system (66 1 in a form which is more convenient for the purpose of computation. In the 
first place, we note that the equations (66 a) and (|66|b) governing the average pressure a can be solved 
by quadrature when the angle 9 satisfies the compatibility conditions for the mixed derivatives of a with 
respect to x and y. This equation takes the form 



'xy 



) cos(2(9) + {9 XX - A9 x 9 y - 9 yy ) sin(20) = 



(67) 



If 9 satisfies (|67|), then the pressure a can be expressed in terms of the velocities u, v, of the angle 9 and 

u(t,x,y) 2 + v(t,x,y) 2 



of the potential V with the form 



1 



v(t, x, y) = -pV(t, x, y) + - sin(20(t, x, y)) + p- 



p / u t (t,x,y)dx 



+ 



(68) 



9 y (t, x, y) sui(29(t, x, y))dx + co(t) 



It should be noted that in order to obtain the compatibility condition (67 1, one has to make use of 



equation (66 e), that is, we must suppose that the flow is irrotational. In order to find a solution of the 
original system ( 66 ) , it is still necessary to find a solution of the overdetermined system consisting of 



equations (67 1, (66 c)-(66 e) for the three unknown functions u,v and 9. In order to reduce the order of 
the PDE (67), we introduce the dependent variables (f> and ijj, which are defined by the equations 



(69) 
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Therefore, we have to solve the following well-determined first-order system of 4 equations in the 4 
unknowns <fi, ip, u, v 



(24? - (</>. 



Vy 



y 1P X ) ~ 2lp 2 ) Uy + 

0, 

0, 



-f 4(f)1P - l^y) U X = 



(70) 



(a) 
(b) 
(<0 

((1) U,y 

Time appears in system ( 70 1 only as a parameter since the system does not involve any derivatives with 
respect to time. Therefore, when solving the system (70 1 we treat the unknowns (j},ip,u,v in (70 1 as 
functions of x and y only, and subsequently we replace the integration constant by arbitrary function of t. 
Since equation (70 a) is not quasilinear, the CSM as presented in section [3] cannot be applied to the entire 
system ( 70 1 . Nevertheless, the CSM can be used for the subsystem A of 3 equations in the 4 unknowns 
consisting of (70 b)-(70 d), which is then underdetermined. Next, we use the degrees of freedom of the 
obtained solution to satisfy equation (70 a) . The subsystem A can be written in the following matrix 
form: 






f 


+ \ 












/ dx 




u 






V 


v J 





10 
1 
1 




(71) 



The dispersion relation (Hi takes the form 



rank 



A 2 
















Ai 


A 2 








A 2 


-Ax 



< 3. 



By solving the dispersion relation ( 72 1 , we obtain the following two types of wave vectors 

A = (l,i), A = (l,-i). 
The Ricmann invariants associated with the wave vectors A and A are 



We seek a real solution of (70 b)-(70 d) of the form 

4> =fi(r, r) +/i(r, r), 



Introducing the notation 



=h{r, r) + / 3 (r,r), 
A = 



r = x — ly. 

i> = h{r, f) 
v = U{r, f) - 

u 4 and 



- fiir, r), 



(72) 

(73) 
(74) 

(75) 



we observe that the matrices r\ a . defined by (43 1 all vanish. Consequently, all trace conditions (49 1 used 
to obtain solutions of the form (75) are identically satisfied. We still have to consider the trace conditions 



(45 1 where the matrices A^, fi — 1,2,3, are given by 
A 1 -- 



-1 

1 








A 1 



10 
1 



A" 



-1 
0-10 
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Conditions ( 45 1 take the form 



(a) 
(b) 
(c) 

Condition (|76"|a) implies that 



df 1 


, dp 


df 1 


dp 


Vi (dp 


, dp 


dp 


dp 


dr 


dr 


dr 


df 


" V dr 


dr 


+ 4 
dr 


df 


dp_ 


dp 


dp 


dp 




y dp_ 


_dp_ 


dp 


dr 


dr 


dr 


dr 


\ dr 


dr 


df 


df 


dp 


dp 


dp 


df 3 


Vi (dP 


dp 


dp 


dp 


dr 


dr 


df 


df 


" V dr 


dr 


+ + 
dr 


df 



(76) 



Tie 



dt]j 
dr 



-Tm 



dr 



(77) 



while conditions (76 b) and (76 c) imply that the solutions for velocities u and v take the particular form 

u = h(r) + h(f), v = i(h(r)-h(f)) , (78) 



where r and f are given by (74 1 and the function h and its complex conjugate are also determined 
using condition (77 1 and the PDE (70 a). This allows us to also determine the functions <fi and tp. For 
computational purposes, it is convenient to introduce velocities of the form (78 1 in equation (66 c) in 
order to determine the angle 9 in the form 



1 / hW(r) 
- arctan i , , . - —. . - 

2 I {hW{r)-hW{?)) 



(79) 



where we denote /iW = d n h/dr n . Next, we substitute (79 1 into equation (67 1 in order to determine the 
function h and its complex conjugate. Proceeding in this way, we find that h and h satisfy the third-order 
ODE which is separable in r and f. It is therefore equivalent to the system 



(a) 



(b) 2 



(80) 



/l(!)(f 



hW(ff 



where fl is a real separation constant. Equation (80 b) can be solved in a way similar way to (80 a). 
Defining 

g(r) = hW(r), (81) 



equation ( 80 a) can be written as 



99 



3/2(g'Y + (n/2)g d = 0. 



(82) 



According to |23| equation ( 82 1 admits a first integral which can be rewritten in the equivalent form 

dg 



VVci-nin(s) 



(83) 



16 



By the change of variable 



exp 



C\ — UJ 



(84) 



equation ( 83 1 can be transformed to 2exp ((20) 1 (w 2 — Ci)) dui — ±ild£, for which the solution, in term 



of the inverse error function erf . is 



= W-2fi 



erf 



e ci/(2H)_ 



-j(c2 ± ng) 

V2ttU 



1 2 



(85) 



where c 2 is a constant of integration. We substitute (84) into (81 ) with uj given by (85 1 and integrate the 
result. We obtain the function h(r) and its complex conjugate in terms of the error functions erfi and 
erf -1 as 

h(r) = — erfi [erf -1 (c 2 + cir)] + c 3 , 

(86) 



h(f) 



2irci 



erfi [erf 1 (c 2 + cir)] + c 3 , 



where the q € C are integration constants and the real separation constant fl is non zero. Since equations 
(86 1 solve the system (80), we have that 9 given by (79) satisfies the compatibility condition (67 1 for a. 
One should note that no derivatives with respect to time appear in the PDE (67 1. Consequently, the 
equation (67 1 is still satisfied if, in equations (78 1 and (79 1, we replace the function h(r) and its complex 
conjugate respectively by 



h(t,r) = - 
h(t,f) = - 



erfi [err 1 (c 2 (t) + Cl (t)r)] + c 3 (t), 
erfi [erf -1 (c 2 (<) + c^f)] + c 3 (t), 



(87) 



which are obtained by substituting the arbitrary complex functions Ci{t) in place of the integration 
constant Ci and the arbitrary real function f2(t) in place of the separation constant tt. Hence, the general 
solution of the system (66) takes the form 



u(t,x,y) = h(t,r) + h(t,r), 
v(t,x,y)=i{h(t,r)-h(t,f)) 



0(t,x,y) 
a(t,x,y) 



arctan i 



hW(t,r) + hW(t,r) ' 
[hW{t,r)-hW(t,f)) 



-pV(t, x,y) + \ sin(20(f , x, y)) + p U(t ' X ' V? ± ^ ^ V? + p / u ,{!..,: !n ,l.r 
9 y (t, x, y) $m(29{t, x, y))dx + a (t), 



where the function h is defined by (87 1, h is its conjugate. Note that the solution (88 1 is a real solution. 



Let us now consider the situation when is identically equal to zero in the system of ODEs ( 80 ) . This 
leads to two possible independent solutions for h(r) depending on whether its second order derivative h rr 



17 



vanishes or not. These solutions are respectively 



(i) h(r) = cir + c 2 , if ^ (2) = 0, 

(**) %) = r^+c 3 , if^ (2) ^0, 



(89) 



r + c 2 

and their complex conjugates, where the Cj are integration constant. Let us consider separately two cases: 



when the function h is given by (89 i) and when h is of the form (89 ii). 



Case i. In this case, the real solution of the original system (66 1 takes the form 



u(t, x, y) =4 {Ke(ci(t))x +Xm( Cl (t))y + Ke(c 2 (t))) 
v(t,x, y) =4 (Im( Cl {t))x - Ke(ci (t))y + Xm(c 2 (*))) 
1 



9(t,x,y) 
a(t,x,y) 



- arctan , , (( 

2 VIm(ci (t)) 



7T 

4 



(90) 



if Xm(ci(t)) ^ 0, 
if Xm(ci{t)) = 0, 

- P V(t, x,y)+p (2 || Cl (0|| 2 + fte(ci(t))) x 2 - 2plm(c 1 (t))xy 

+ p (2 || Cl (<)|| 2 - ftefaft))) y 2 + 2p (2fte (ci(i)pa(i)) + ^e(c 2 (t))) x 

- 2p (21m (ci(t)c 2 (t)) +Xm(c 2 (t))) y + a (t), 

where the real function tro(i) and the complex functions C\(t), c 2 (t) are arbitrary functions of time and 
Ci(t), c 2 (t) their respective derivatives with respect to t. It should be noted that if the potential V, the 
real function ao(t) and the complex functions ci(t), c 2 (t), are all bounded and are their derivatives C\(t), 



c 2 (t) are also bounded, then the solution (83) is bounded. Moreover, if we take these functions to be of 
the form 



-s t 



Cl (t) = a^e"^* + ib je - q i\ j = 1, 2, 



(91) 



where < qj,Sj g M, a,j,bj G BL j 
interesting situation occurs when the arbitrary functions Cj (t) are constant in solution ( 83 1 . In this case 



0,1,2, then we obtain a damped solution. Another physically 



we obtain a stationary solution for velocities u and v. This allows us to draw the shape of those extrusion 
dies which are admissible by requiring that the tool walls coincide with the lines of flow generated by the 
velocities u and v. This is necessary since we suppose that the flow is incompressible. From the practical 
point of view, it is convenient to press the plastic material through the die rectilinearly with constant 
speed. If we feed the die in this way, the plasticity region limits are curves defined by the ODE 



dy_ _ Vq - v(x,y) 
dx Uq — u(x, y) ' 



(92) 



where Uq and Vq are constants representing the feeding velocity (or extraction velocity) of the die along 



the cc-axis and y-axis respectively. Equation ( 92 1 is a consequence of the hypothesis that the flow is 



incompressible and the mass is conserved. This condition reduces to the boundary conditions described 
in [S] when we require that the limits of the plasticity region correspond to the slip lines (which correspond 
to the characteristic curves of the original system ( 66 1 ) , that is, when we require that dy/dx = tan# or 
dy/dx = —cot 6. Here, we use the weakened condition (92 1 because no constraints are imposed on 
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Figure 1 : Extrusion die corresponding to the solution ( 83 1 



the flow lines which are in contact with the tool walls. This is so because, for a given solution and 
given parameters, we choose the walls of the extrusion die to lie along the flow lines. For the purpose 
of illustrating the applicability of the method, we have drawn in figure ^ the shape of a tool and the 
flow of matter inside the extrusion die for the following parameters : lZe(ci(t)) = 1, Tm(c\{t)) = 0, 
lZe(c 2 (t)) = 0, Xm{c 2 {t)) = 0. The feeding velocity of the tool is Uo — 5.95, Vq — and the tool expels 
the matter at a velocity of U\ = 24.05, V\ = 0. The plasticity region at the opening is bounded by curve 
C\ and at the exit by the curve C 2 . This extrusion die can thin a plate or rod of ideal plastic material. 



Case ii. If h is defined by (82 ii), then the corresponding solution of system ( |66[ ) takes the form 

Tie (ci (t)pa (t) ) + Ke{ Cl {t))x + Xm( Cl {y))y 



u{t,x,y) =Ke(c a (t)) 



v{t,x,y) =He(c 3 {t)) + 2 



x 2 + y 2 + 2 {1Ze(c 2 (t))x+lm(c 2 (t))y) + \\c 2 {t)\Y 

-1m (a(t)c 2 (t)) - Xm(ci(£))x + Kejcxjy^y 
x 2 + y 2 + 2 {Ke(c 2 (t))x +Xm(c 2 (t))y) + \\c 2 (t)f 

lm( Cl (t))(x 2 - y 2 ) + 2Ue(c 1 (t))xy 



6{t, x, y) = — — arctan 

+ (-lm{ Cl {t)c 2 (t))+ne(cx{t))ne{c 2 (t)) -lm{c x {t))ne{c 2 (t)))x 
+ (Ke{c 1 {t)c 2 {t))+'Re{c 1 {t))'Re{c 2 {t)) +lm{c 1 {t))'Re{c 2 {t)))y 

+ He (ci (t)c 2 (t)) Ue(c 2 (t)) + 1m (a (t)c 2 (t))lm(c 2 (t)) 
7le(ci(t))(a; 2 - y 2 ) + 2Xm(c 1 (t))xy 

+ ( - TZe( Cl (t)c 2 (t)) +Ke(c 1 (t))Ke(c 2 (t)) +lm{ Cl {t))Tm{c 2 (t)))x 
+ (lm( Cl (t)c 2 (t)) - ne{ Cl {t))Tm{c 2 {t)) +Tm( Cl (t))Ke(c 2 {t)))y 

+ Ke { Cl (t)c 2 (t))1le(c 2 {t)) +Xm (c 1 {t)c 2 {t))Xm{c 2 {t)) 



(93) 



19 



Figure 2: Extrusion die corresponding to the solution (85 1 



where the mean pressure a is given by (68 1, where we substitute the values of functions u,v,9 given by 
(85 1. The complex functions Cj(i), i — 1,2,3, and the real function ao(t) which appear in this solution 
are arbitrary. For any time to, solution ( |68[ ), ( |93[ ) has a singularity at point (xo,yo) which satisfies the 
equation 

xl + yl+2 {Ke(c 2 {to))x + lm(c 2 {t ))y a ) + \c 2 {t )\ 2 = 0. 
This singularity is stationary at the origin if the function c 2 (t) is constant. Otherwise, its position varies 



with time. If the functions do and Cj, i = 1, 2, 3, are of the form (91 1 defined in a region of the a;y-plane 
on a time interval \Tq,T) where the gradient catastrophe does not occur, then the solution is bounded 
and damped. In figure [2] we have drawn the shape of an extrusion die corresponding to the solution ( 93 1 



for the following choice of parameters: IZe (ci(t)) = 0, Ira (ci(£)) = 0, IZe (c 2 (t)) = 0, Ira {c 2 (t)) = —0.5, 
TZe(cs(t)) = —0.5 and Im{c^{t)) = 0. The feeding velocity has component Uq — 0.2, Vq — 0.2, and the 
extraction of material is performed at the velocity U\ — 0.2, V\ — —0.2. This type of tool can be used to 
bend a rod by extrusion without having to fold it. 



Finally, we should emphasize that the flow changes considerably when the parameters are varied and 
we have the freedom to choose the walls of the tool among the flow lines for certain fixed parameters 
Cj(i). Moreover, the velocity and the orientation of the feeding (extraction) can vary somewhat for a 
given shape of the tool. Consequently, many types of extrusion dies can be drawn. 



5 Final remarks. 

The generalized method of characteristics was originally devised for solving first-order quasilinear hy- 
perbolic systems. The proposed technique allows us to extend the applicability of this approach not 
only to hyperbolic systems but also to encompass elliptic and mixed (parabolic) type systems. A vari- 
ant of the conditional symmetry method for obtaining multimode solutions has been proposed for these 
types of systems. We have demonstrated the usefulness of this approach through several examples of 
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the ideal plasticity in (2 + 1) dimensions. New classes of real solutions have been constructed in closed 
form, some of them bounded. Some of the obtained solutions have a stationary flow for an appropriate 
choice of parameters Cj(i). For these solutions, we have drawn extrusion dies and the vector fields which 
define the flow inside a region where the gradient catastrophe does not occur. The proposed approach 
for constructing multimode solutions can be used in several potential applications arising from systems 
describing nonlinear phenomena in physics or in theoretical chemistry and biology. It should be noted 
that in the multidimensional case for many physical models there are very few known examples of multi- 
mode solutions. The approach proposed in this paper offers a new and promising way to construct and 
investigate such types of solutions. This makes our approach attractive since it can widen the potential 
range of applications leading to more diverse types of solutions. This task will be undertaken in a future 
work. 
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